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AN ACTIVE CURVE APPROACH FOR TOMOGRAPHIC RECONSTRUCTION OF 
BINARY RADIALLY SYMMETRIC OBJECTS 

I. Abraham^ R. Abraham and M. Bergounioux^ 

Abstract. This paper deals with a method of tomographic reconstruction of radially symmetric 
objects from a single radiograph, in order to study the behavior of shocked material. The usual 
tomographic reconstruction algorithms such as generalized inverse or filtered back-projection can- 
not be applied here because data are very noisy and the inverse problem associated to single view 
tomographic reconstruction is highly unstable. In order to improve the reconstruction, we propose 
here to add some a priori assumptions on the looked after object. One of these assumptions is that 
the object is binary and consequently, the object may be described by the curves that separate the 
two materials. We present a model that lives in BV space and leads to a non local Hamilton- Jacobi 
equation, via a level set strategy. Numerical experiments are performed (using level sets methods) 
on synthetic objects. 
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1. Introduction 

Medical scanner is the most used application of tomographic reconstruction. It allows to explore 
the interior of a human body. In the same way, industrial tomography explores the interior of an 
object and is often used for non-destructive testing. 

We are interested here in a very specific application of tomographic reconstruction for a physical 
experiment described later. The goal of this experiment is to study the behavior of a material under 
a shock. We obtain during the deformation of the object an X-ray radiography by high speed 
image capture. We suppose this object is radially symmetric, so that one radiograph is enough to 
reconstruct the 3D object. 
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Figure 1 . Experimental setup 



Several authors have proposed techniques (standard in medical tomography) for tomographic 
reconstruction when enough projections (from different points of view) are available: this allows 
to get an analytic formula for the solution (see for instance [9] or [6]). These methods cannot 
be used directly when only a few number of projections is known. Some alternative methods 
have been proposed in order to partially reconstruct the densities (see for instance [5]). We are 
interesting here in single view tomographic reconstruction for radially symmetric object (see for 
instance [8] for a more complete presentation of the subject). As any tomographic reconstruction, 
this problem leads to an ill-posed inverse problem. As we only have one radiograph, data are not 
very redundant and the ill-posed character is even more accurate. 

We present here a tomographic method adapted to this specific problem, originally developed 
in [1], and based on a curve evolution approach. The main idea is to add some a priori knowledge 
on the object we are studying in order to improve the reconstruction. The object may then be 
described by a small set of characters (in this case, they will be curves) which are estimated by the 
minimization of an energy functional. This work is very close to another work by Feng and al [7]. 
The main difference is the purpose of the work: whereas they are seeking recovering textures, we 
are looking for accurate edges. It is also close to the results of Bruandet and al [3]. However, the 
present work handles very noisy data and highly unstable inverse problems, and shows how this 
method is powerful despite these perturbations. Further, we take here into account the effects of 
blur (which may be non-linear) and try to deconvolve the image during the reconstruction. 

Let us mention at this point that our framework if completely different from the usual tomo- 
graphic point of view, and usual techniques (such as filtered back-projection) are not adapted to 
our case. Indeed, usually, as the X-rays are supposed to be parallel (this is also the case here), 
the "horizontal" slices of the object are supposed to be independent and studied separately. Usual 
regularization techniques deal with one slice and regularize this particular slice. Here, because of 
the radial symmetry, the slices are composed of concentric annulus and do not need any regular- 
ization. The goal of this work is to add some consistency between the slices in order to improve 
the reconstruction. 

The paper is organized as follows. First we present the physical experiment whose data are 
extracted and explain what are the motivations of the work. Next, we introduce the projection 
operator. In Section 4, we present a continuous model with the suitable functional framework and 
prove existence result. Section 5 is devoted to formal computation of the energy derivative in order 
to state some optimality conditions. In Section 6, a front propagation point of view is adopted and 
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the level set method leads to a non local Hamilton-Jacobi equation. In the last section, we present 
some numerical results and give hints for numerical schemes improvement. 

2. Experiment 

This work is part of some physical experiments whose goal is the study of the behavior of 
shocked material. The present experiment consists in making a hull of well known material im- 
plode using surrounding explosives. The whole initial physical setup (the hull, the explosives 
...) are radially symmetric. A reasonable assumption is to suppose that during the implosion, 
everything remains radially symmetric. 

Physicists are looking for the shape of the interior at some fixed time of interest. At that time, 
the interior may be composed of several holes which also may be very irregular. Figure |3l is a 
synthetic object that contains all the standard difficulties that may appear. These difficulties are 
characterized by: 

• Several disconnected holes. 

• A small hole located on the symmetry axis (which is the area where the details are difficult 
to recover). 

• Smaller and smaller details on the boundary of the top hole in order to determine a lower 
bound detection. 

To achieve this goal, a X-rays radiograph is obtained. In order to extract the desired informa- 
tions, a tomographic reconstruction must be performed. Let us note here that, as the object is 
radially symmetric, a single radiography is enough to compute the reconstruction. 

A radiography measures the attenuation of X-rays through the object. A point on the radio- 
graphy will be determined by its coordinates {u, v) in a Cartesian coordinates system where the 
f -axis will be the projection of the symmetry axis. If Iq is the intensity of the incident X-rays flux, 
the measured flux / at a point (n, v) is given by 

/ = Joe" ^ ^^i■'^0'^)d<^ 

where the integral operates along the ray that reaches the point (u, v) of the detector, d£ is the in- 
finitesimal element of length along the ray and fj, is the local attenuation coefficient. For simplicity, 
we will consider that this coefficient is proportional to the material density. To deal with linear 
operators, we take the Neperian logarithm of this attenuation and will call the transformation 




the projection operator. 

Through the rest of the paper, in order to simplify the expression of the projection operator, we 
will suppose that the X-ray source is far enough away from the object so that we may consider 
that the rays are parallel, and orthogonal to the symmetry axis. As a consequence, the horizontal 
slices of the object may be considered separately to perform the projection. 

As the studied object is radially symmetric, we will work in a system of cylinder coordinates 
(r, 6, z) where the z-axis is the symmetry axis. The object is then described by the density at the 
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point (r, 9, z), which is given by a function / which depends only on (r, z) by symmetry. In the 
text, the notation / will always refer to the density of the object. A typical function / is given in 
Figure|2l It represents an object composed of concentric shells of homogeneous materials (called 
the "exterior" in what follows) surrounding a ball (called the "interior") of another homogeneous 
material that contains some empty holes. This figure may be viewed as a slice of the object by a 
plane that contains the symmetry axis. To recover the 3D-object, it suffices to perform a rotation 
of this image around the z axis. For instance, the two round white holes in the center are in fact 
the slice of a torus. As the looked-after characteristic of the object is the shape of the holes, we 
will focus only on the interior of the object (see Figure |3l). We here handle only binary objects 
composed of one homogeneous material (in black) and some holes (in white). 




Figure 2. Slice of a typical binary radially symmetric object by a plane that 
contains the symmetry axis (the z-axis). 




Figure 3. Zoom on the interior of the object of Figure |2l The homogeneous 
material is in black whereas the holes are in white. 
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3. A VARIATIONAL APPROACH 

3.1. The projection operator 

We first explicit the projection operator and its adjoint. 

Proposition 3.1. In the case of a radially symmetric object, the projection operator, denoted by 
H, is given, for every function f G L°°(R+ x M) with compact support, by 

r+oo 

V(«, t-) G R X M Hf{u, v) = 2 f{r, v) _ _ dr. (3.1) 

J\u\ vr^ — ii 

Proof - Consider a 3D-object which is described by a function f{x, y, z) (Cartesian coordinates 
system). The projection operator H is 



Hf{u,v)= / f{x,u,v)dx. 
Jr 



In the case of radially symmetric object, we parametrize the object by a function f{r,z) with 
cylinder coordinates. Therefore 



f{x,y,z) = fiV^^+y^,z) . 

Then, we have, for u > 

r ^ r r+oo 

Hf{u,v)= / f{x,u,v)dx= / /(V x'^ + u^jv) dx = 2 / /(yx^^i^, 'u) cJa; . 
Jr Jr Jo 

We perform the following change of variable 

r = x >0 X = \J — xP' ^ r > u , 



to get 



For u < 0, we have 



r+oo 

Hfiu,v) = 2 / fir,v) ^ dr 

Ju V — u 



Hf{u,v) = 2 r f{r,v) , i''' =dr 

J —oo 



r\ 

2 — 7,2 



\fr^ — u 



Using the change of variable u:to — u and the fact that r f{r, v) is even, by symmetry, we get 

Hf{u,v) = 2 f{r,v)^^=dr. 

J\u\ \/r'^ — u^ 



\u\ 

□ 
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Remark 3.1. Operator H may be defined by density on measurable fiinctions f such that all the 
partial applications /(•, z) belong to L^(M+). Then, all thefimctions Hf{-,v) belong to the space 
BMO(K) of bounded mean oscillations functions : 

SMO(R) = {/ : R ^ R I sup ( 4 / \f{x) - dx] <+oo}, 

R>0 \^ J\x-y\<R J 

I ry+R 

where fR{y) = / f{x) dx. For more details, one can refer to [13]. 

2-fl Jy-R 

In the sequel, we will need to handle functions / that are defined on R^ (instead of on R_)_ x R). 
We thus define the operator H for function / G L°°(R^) with compact support by 



r+oc J, 
Hf{u,v) = 2 f{sgn{u)r,v) dr 
J\u\ vr^ — w^ 

although this has no more physical meaning. Here, the function sgn is defined by 

sgn{x) 



1 ifx>0, 
-1 ifa;<0. 



We shall also need the back-projection that is the adjoint operator H* of H; it can be computed 
in a similar way. 

Proposition 3.2. The adjoint operator ( in L?) H* of the projection operator is given, for every 
function g G L°°(R^) with compact support by : 

f\r\ 

Vr G R , Vz G R, H*g{r, z) = 2 g{sgn{r)u, z) ' du. (3.2) 

Jo Vr — V? 

Proof - The adjoint operator H* of H is the unique operator such that, for every / and g in 
L°°(R^) with compact support, 

/+0O r+oo r+oo r+oo 

/ Hf{u,v)g{u,v)dvdu= / / f{r,z)H*g{r,z)drdz. 
-oo J —oo J —oo J —oo 
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Using (I3.1I) and Fubini's theorem, we get 



■oo ^+00 

Hf{u, v) g{u, v) dv du 
00 J —00 

'u=0 /•v=+oo f+00 



u=—oo Jv=—oo Jr=\u\ \/r^ U 



r 

f{-r, v)^===g{u, v) dv du dr 



u=+oo ry=+oo r+00 ^ 



+ 2 / / f{r,v)—===g{u,v)dvdudr 

J u=0 J v=—oo J r=u v r U 



v=+oo fr=+oo fu=0 ^ 



v=—oo Jr=0 Ju=—r \fr^ 



/(— r, v)g{u, v) du dr dv 



D=+oo />r=+oo cu=r ^ 



+ 2 / / -^==f{r,v)g{u,v)dudrdv 

'v=—oo Jr=0 Ju=0 \/r — U 



ti=+oo fr=— 00 fu=0 ^ 



-f{r, v)g{u, v) du (—dr) dv 



v=—oo J r=0 J u=r \fr^— U^ 

/•v=+OD /•r=+oo i'U=r ^ 

+ 2 / / f{r,v)g{u,v)dudrdv 

Jv=-oo Jr=o Ju=o yr-^ — 

2 I / / —===f{r,v)g{u,v)dudrdv 

v=—oo Jr=— 00 J u=0 \r U 

2 / / / , ^ f{r,v)g{sgn{r)u,v)dudrdv 



v=—oo Jr=—oo J u=0 \/ r'^ U?' 

So we obtain the following expression for the back projection : 

|r| 



ii*g{r,z)=2 [ 
Jo 



r\ 

g{sgn{r)u, z) du 



_ y2 



3.2. Toward a continuous model 



□ 



Thanks to the symmetry, this operator characterizes the Radon transform of the object and so 
is invertible; one radiograph is enough to reconstruct the object. The inverse operator is given, for 
an almost everywhere differentiable function g with compact support, and for every r > 0, by 



H-'g{r,z)= / o x-^' ^ ' dx. 

Jr yx'^ — r'^ 

Because of the derivative term, the operator H^^ is not continuous. Consequently, a small vari- 
ation on the measure g leads to significant errors on the reconstruction. As our radiographs are 
strongly perturbed, applying H^^ to our data leads to a poor reconstruction. Due to the experi- 
mental setup they are also two main perturbations: 

• A blur, due to the detector response and the X-ray source spot size. 
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• A noise. 

Others perturbations such as scattered field, motion blur... also exist but are neglected in this study. 
We denote by F the effect of blurs. We will consider the following simplified case where F is 
supposed to be linear 

F(k) = N*k (3.3) 

where * is the usual convolution operation, k is the projected image and is a positive symmetric 
kernel. 

Remark 3.2. A more realistic case stands when the convolution operates on the intensity : then 
F is of the form 

F{k) = ^ In (^e-^'' * iv) 

where u is the multiplicative coefficient between the density and the attenuation coefficient. Some 
specific experiments have been carried out to measure the blur effect. Consequently, we will 
suppose that, in both cases, the kernel N is known. The linear blur is not realistic but is treated 
here to make the computations simpler for the presentation. 

The noise is supposed for simplicity to be an additive Gaussian white noise of mean 0, denoted 
by e. Consequently, the projection of the object / will be 

g = F{Hf) + e. 

The comparison between the theoretical projection F{Hf) and the perturbed one is shown on 
Figure |4] The reconstruction using the inverse operator H^^ applied to g is given by Figure |5] 
The purpose of the experiment is to separate the material from the empty holes and consequently 
to precisely determine the frontier between the two areas, which is difficult to perform on the 
reconstruction of Figured 




Figure 4. Left-hand side: theoretical projection F{Hf) of the object of Fig- 
ure |5] Right-hand side: real projection of the same object with realistic noise and 
blur. 
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Figure 5. Comparison between the real object on the left-hand side and the 
reconstruction computed with H^^ applied to the real projection on the right-hand 
side. 

It is clear from Figure |5lthat the use of the inverse operator is not suitable. In order to improve 
the reconstruction, we must add some a priori knowledge on the object to be reconstructed. Indeed 
the object that we reconstruct must satisfy some physical property. 

We chose to stress on two points: 

• The center of the object is composed of one homogeneous known material's density with 
some holes inside. 

• There cannot be any material inside a hole. 

In a previous work [10], J.M. Lagrange reconstructed the exterior of the object. In this recon- 
struction, the density of the material at the center of the object is known, only the holes are not 
reconstructed. In other words, we can reconstruct an object without holes and we can compute (as 
H and F are known) the theoretical projection of this reconstruction. We then act as if the blurred 
projection was linear and subtract the projection of the non-holed object to the data. In what 
follows, we will call experimental data this subtracted image which corresponds to the blurred 
projection of a "fictive" object of density with some holes of known "density" A > 0. Conse- 
quently, the space of admitted objects will be the set of functions / that take values in {0, A}. This 
space of functions will be denoted JP" in the sequel. 

The second hypothesis is more difficult to take into account. We chose in this work to tackle 
the problem via an energy minimization method where the energy functional is composed of two 
terms: the first one is a matching term , the second one is a penalization term which tries to handle 
the second assumption. The matching term will be a L^-norm which is justified by the Gaussian 
white noise. 

Remark 3.3. In the case where F is not linear, the exact method to remove the exterior is to 
operate the blur function on the addition of the known exterior and the center. For the sake 
of simplicity, we will not use this method and will consider that the errors are negligible when 
subtracting the projections. 

Let us first describe more precisely the set JT. This time, the functions f ^ T will be defined 
on M^, with values still in {0, A}, with compact support. Therefore, such a function / will be 
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characterized by the knowledge of the curves that limit the two areas where / is equal to A and 
to 0. Indeed, as the support of the function / is bounded, these curves are disjoint Jordan curves 
and the density of the inside is A whereas the density of the outside is 0. Consequently, the energy 
that we will consider will be a function of 7^ where 7^ is a set of disjoint Jordan curves. For 
mathematical reasons, we must add an extra-assumption : the curves 7^ are so that the normal 
vector of the curves is well-defined (as an orthogonal vector to the tangent one). 

In this continuous framework, the matching term is just the usual L^-norm between Hf and the 
data g (where H is given by (13.11) ). So, the first term is 

E,i^f) = \\FiHf)-g\\l 

For the penalization term, we choose 

^2(7/) = ^(7/) 

where ^(7/) denotes the length of the curves jj. Let us remark that this penalization term may be 
also viewed as the total variation (up to a multiplicative constant) of the function / because of the 
binarity. Eventually, the total energy functional is 

E{jf) = \\F{Hf)-g\\l + aei^f) (3.4) 

which is an adaptation of the well-known Mumford-Shah energy functional introduced in [11]. 
The "optimal" value of a may depend on the data. 

3.3. A continuous model in BV space 

The previous analysis gives the mains ideas for the modelization. Now, we make it precise 
using an appropriate functional framework. Let 17 be a bounded open subset M? with Lipschitz 
boundary. We shall consider bounded variation functions. Recall that the space of such functions 
is 

BV{n) = {u£ L^{n) I j(u) < +00} 

where 

J{u) = snj> !^j^u{x) div ^{x)dx \ ^ eCUn) , ||e||oo<l|, (3.5) 

where C^(il) denotes the space of functions with compact support in Q. The space BV{^) 
endowed with the norm 

= IKIIli + Jiu) , 

is a Banach space. 

If u G BV{^1) its derivative in V^Q) (distributions) is a bounded Radon measure denoted Du 
and J(n) is the total variation of \Du\ on 17. Let us recall useful properties of BV-functions ( [2]): 

Proposition 3.3. Let Q be an open subset with Lipschitz boundary. 
L If u G BV{il), we get the following decomposition for Du : 



Du = Vudx + D^u , 
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where Vudx is the absolutely continuous part of Du with respect of the Lebesgue measure and 
D^u is the singular part. 

2. The map u ^ J {u) from BV{Q.) to R+ is lower semi-continuous (Isc)for the L^(ri) topology. 

3. BV(^l) C with compact embedding. 

4. BV{p,) C L^(J7) with compact embedding. 

We precise hereafter an important continuity property of the projection operator H. 

Proposition 3.4. The projection operator H is continuous from L^+*(f2) to LP {^) for every p G 
[1, +oo] and s > 0. 

Proof - It is a direct consequence of Holder inequality. Let be s > and / G L^'^^{Q,). Its 
support is included in Q which is included in some [— M, +M] x [— M, +M] where M > and 
only depends on il. It is clear that Hf is defined everywhere on Cl and, for every u > 



\Hf{u,v) 



+ CXD 



f{r,v) 



^Jr^ — u' 



; dr 



M 



fir,v) 



y/r"^ — u' 



; dr 



< 2 



M 



\f{r,v)f+'dr 



1 

2+s 



Ju 



(r + u) 2 (r — u) 



dr 



where q = 1 + 



Therefore 



1 + s 

\Hfiu,v)\<2\\f\\L2+s 



M 



r2[r — u)2 



-dr 



<2M-2\\f\\L2+s 



M 



{r — u) 2 dr 



The computations in the case n < are similare and lead to the same inequality (with some 
additional absolute values). As 



1-5 



we get 



2 2(1 + s) 



\Hf{u,v)\ < 2Mi\\f\\L2+s ^^l±ll [M-u]^ < Ci^},s)\\f\\L2+s 



(3.6) 



here and in the sequel C(ri, s) denotes a generic constant depending on s and So 

\\Hf\\^<c{n,s)\\f\\L2+s. 

As Q, is bounded, this yields 

V/ G L^+%n), Vp G [l,+oo] \\Hf\\L.^n) < C{n,s)\\f\U.+s . 



(3.7) 

□ 
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As we consider the length of curves, the most suitable functional space to set a variational 
formulation of the reconstruction problem is BV{^). Therefore, we consider the following mini- 
mization problem 

r mm\\FHf-gg + aJif) 

{V) I feBV{n) 

|/(x)| = 1 a.e. on Q, 

Here 

• II • II2 stands for the L^(r2)- norm, g G L^(i7) and a > 0. 

• The operator F is given by (I3.3t . Without loss of generality, we may assume (for simplic- 
ity) that F = I. 

• At last, "|/(x)| = 1 a.e. sur 0", is the binarity constraint. We have mentionned that 

the image fo takes its values in {0, A} where A > 0. With the change of variable 
2 

/ = — — /o + 1, we may assume that the image values belong to {—1, 1}. 
A 

Remark 3.4. A similar problem has been studied in [4] with smoother projection operator and 
convex constraints. This is not our case. The pointwise constraint " |/(x)| = 1 a.e. on Vl" is 
a very hard constraint. The constraint set is not convex and its interior is empty for most usual 
topologies. 

Now we may give the main result of this section : 
Theorem 3.1. Problem (V) admits at least a solution. 

Proof - Let ipn G BV{Q) be a minimizing sequence. It satisfies Hy'nlloo = 1 ; so 

Vp G [l,+oo[, Vn G N, W^PtiWlp < ■ (3.8) 

Therefore the sequence is L^(i7)- bounded. As J{ipn) is bounded as well, the sequence is 
bounded in BV{Q,). Thus it converges (extracting a subsequence) to some ip G BV{^1.) for the 
weak-star topology. 

Estimate (13.81) implies the weak convergence of {ipn) to ip in L'^~^^{^) for every s > 0. Thanks 
to the H continuity property of proposition 13.41 we assert that Hipn weakly converges to Htp in 
L^{n). We get 

\\Hip - g\\l2 < liminf \\Hipn - gWl^ , (3.9) 
with the lower semi-continuity of the norm. 

Moreover BV{Q) is compactly embedded in L^(i7). This yields that strongly converges to 
if in L^(0). As J is Isc with respect to L^{Q)- topology, we get 

J(V5) < liminf J((^„) , (3.10) 

n— *oo 

Finally 

inf(P) = lim \\Hipn - gWli + aJ{ipn) 

> liminf \\Hipn - §11% + aJ{ipn) > \\Hip - gW^z + aJ{ip) . 

n^oo 
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As the pointwise constraint is obviously satisfied, ipisa solution to {V). □ 

4. Computation of the energy derivative 

Now we look for optimality conditions. Unfortunately we cannot compute easily the derivative 
of the energy in the BV{Q) framework. Indeed we need regular curves and we do not know if the 
BV{^) minimizer provides a curve with the required regularity. Moreover, the set of constraints 
is not convex and it is not easy to compute the Gateaux- derivative (no admissible test functions). 

So we have few hope to get classical optimality conditions and we rather compute minimizing 
sequences. We focus on particular ones that are given via the gradient descent method inspired 
by [11]. Formally, we look for a family of curves (7^)^ > such that 

f <-> <- » 

so that £'(7t) decreases as t — > +00. Let us compute the energy variation when we operate a small 
deformation on the curves 7. In other word, we will compute the Gateau derivative of the energy 
for a small deformation 6^: 

dj^'^ ' t^o t 

We will first focus on local deformations ^7. Let (rg, zq) be a point P of 7. We consider a 
local reference system which center is P and axis are given by the tangent and normal vectors at 
P and we denote by rj) the new generic coordinates in this reference system. With an abuse of 
notation, we still denote /(^, rj) = f{r, z). We apply the implicit functions theorem to parametrize 
our curve: there exist a neighborhood U of P and a function h such that, for every rj) G U, 

G 7 V = HO- 

Eventually, we get a neighborhood U of P, a neighborhood / of ^0 and a function h such that 

The local parametrization is oriented along the outward normal n to the curve 7 at point P (see 
figure|6ll. More precisely, we define the local coordinate system (r, n) where r is the usual tangent 
vector, n is the direct orthonormal vector; we set the curve orientation so that n is the outward 
normal. The function / if then defined on U by 



A ifr?</i(0 
ifr?>/i(0 



This parametrization is described on figure |6l We then consider a local (limited to U) deformation 
^7 along the normal vector. This is equivalent to handling a function 6h whose support is 
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included in I. The new curve obtained after the deformation td'j is then parametrized by 

V -- 



h{0+tSh{0 for(e,r/)GC/ 
7 otherwise 



This defines a new function fi 



ft{^,v) = { A if{^,7])eun{v<h{C) + t6h{c)} (4.11) 
if (^,r,)eun{v> HO + t6h{c)} 

We will also set 5ft = ft — f- This deformation is described on figure|6l 




Figure 6. Description of a local deformation of the initial curve 7. P is the 
current point, U is the neighborhood of P in which the deformation is restricted 
to and 7 + is the new curve after deformation. The interior of the curve is the 
set where / = A 



The Gateau derivative for the energy E2 has already been computed in [1 1] and is 

^(7)<57 = - J^cmv{j){C,h{0)5h{0d^ 

where curv denotes the curvature of the curve and 5h is the parametrization of ^7. 
It remains to compute the derivative for the matching term. First we estimate 6ft : a simple com- 
putation shows that 

dft{^,v} - I A if HO <v< HO + t5H0 



and 



Sn(^v)-l ° ifv<HO + t6HOorviO>HO inease<5/.<0 

oit[i;,V) - I -A if HO >v> HO + t^HO 
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Now we compute -E'i(7t) — -^1(7) where 7 (resp. 7t) is the curve associated to the function / 
(resp. ft): 

= [ [ {{g- FHftf -{g- FHff) {u, v) du dv 

= 11 {{g-FHf-FH5ftf-{g-FHff){u,v)dudv 

JR JR 

= -2/ [ {g-FHf){u,v)FH6ftiu,v)dudv+ [ [ {FH5ftf{u,v) dudv 
Jr. Jr Jr Jr 

. ' 

=o{t) 

= ~2{{g-FHf),FH6ft)L2+o{t) 
= -2 {H*F*{g - FHf),5ft)L2 + o{t). 

To simplify the notations, we denote by Af := {H*F*g — H*F*FHf) so that we need to 
compute 

As 6 ft is zero out of the neighbourhood U, we have 

{Af,6ft)L^= [ Af{^,v)Sft{^,v)d^dv. 
Ju 

In the case 6h > 0, we have, 

{Af,Sft)L^ = X / Afi^,r,)d^d7^. 



As the function Af is continuous (and thus bounded on U), we may pass to the limit by domi- 
nated convergence and get 



Ymi^{Af,6ft)L2 = xl^Af{^,h{0)6h{0dC. 



t- 

In the case 6h < 0, we have 

l^el J r,=h(0+tSh(i) 



{Af, 5ft)L^ = / (- A) / / Af{i, v)d^ dn 

J J eel Jv=h(£)+tSh(f) 



and we obtain the same limit as in the nonnegative case. 
Finally, the energy derivative is 



1^ (7/) • S^f = -2\jAf{^, h{0) ShiO dC-»l curv(7/) h{0)Sh{0 d^ 
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If we set /3 = — , we get 



— {jf) . 5jf = -2 (XAf + /3curv(7^)(C, h{0)) 5h{0 . (4.12) 
As 6h =< d-yfjU > formula (I4.12I) may be written 

^(7/)-'57/ = -2 j (A^/ + /3curv(7y)(s)) <%,n> ds (4.13) 



where n denotes the outward pointing normal unit vector of the curve 7, < • , • > denotes the usual 
scalar product in and c{s) is a positive coefficient that depends on the curvilinear abscissa s. 

The latter expression is linear and continuous in ^7, this formula is also true for a non-local 
deformation (which can be achieved by summing local deformations). 



5. Front propagation and level set method 

The goal of the present section is to consider a family of curves (7t)f>o that will converge 
toward a local minimum of the functional energy. From equation (I4.13t . it is clear that if the 
curves (7^) evolve according to the differential equation 

1^ = (A^/ + /3curv(7/))n, (5.14) 

the total energy will decrease. 

To implement a numerical scheme that discretizes equation ( I5.14t . it is easier to use a level set 
method (see [12] for a complete exposition of the level set method). Indeed, equation (I5.14t may 
present some instabilities, in particular when two curves collide during the evolution or when a 
curve must disappear. All these evolutions are handled easily via the level set method. 

The level set method consists in viewing the curves 7 as the 0-level set of a smooth real function 
(j) defined on M^. The function / that we are seeking is then just given by the formula 

f{x) = Al0(a.)>o. 

We must then write an evolution PDF for the functions 0^ = 4){t,-) that corresponds to the curves 
7^. Let x{t) be a point of the curve 7^ and let us follow that point during the evolution. We know 
that this point evolves according to equation 15.141 

x'{t) = {XAf + (3 curv(7/)) {x{t))n. 

We can re-write this equation in terms of the function (p recognizing that 
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where V stands for the gradient of (p with respect to | • | denotes the euclidean norm. The 
evolution equation becomes 

-'(t) = (a^(A1,(,,)>o) +/3div (^)) 

Then, as the point x{t) remains on the curve 7^, it satisfies cl)t{x{t)) = (j){t,x{t)) = 0. By 
differentiating this expression, we obtain 

^ + (V^,x'(t)> = 
which leads to the following evolution equation for (/>: 

^ + I V.0I (a^ (A1,(,,)>o) + /3div (^) ) = , 

that is 

^ = \V4>\ {>?H*F*FH (l<^(i,)>o) - /3div - \H*F*g^ . (5.15) 

The above equation is an Hamilton- Jacobi equation which involves a non local term (through 
H and F). Such equations are difficult to handle especially when it is not monotone (which is the 
case here). In particular, existence and/or uniqueness of solutions (even in the viscosity sense) are 
not clear. The approximation process is not easy as well and the numerical realization remains a 
challenge though this equation is a scalar equation which is easier to discretize than the vectorial 
one. Here we used the discrete scheme described in [12] to get numerical results. 

6. Results and discussion 

6.1. Explicit scheme 

We briefly present the numerical scheme. We used an explicit scheme in time and the spatial 
discretization has been performed following [12] . We set 

G = {FH)*FH and g* = H*F*g , 

so that the equation (15.151) is 

We set tn = nAt, <I>" = (j){tn, ■), X = {xi,yj)(^ij)^i with Xi = iAx et yj = jAy. The explicit 
Euler scheme gives : 

= $"(X) + At |V^>"|(X) (A2g(l$n>o) -/?curv($") - \g*itn,X)) . 
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The curvature term is computed as 



where stands for the partial derivative with respect to x. The discrete approximation of the 
gradient is standard: 

are Dy^ defined in the same way. A usual approximation for |V$| is given by : 

|V$|(X)~ [max(D+$,0)2 + max(D+$,0)2 +min(i:i-$,0)2 + min(D-$,0)2]^/^ . 
The non local term ^(l$n>o) is exactly computed. 

6.2. Numerical results 

The previous scheme has been implemented on a 3.6 GHz PC. A classical reinitialization pro- 
cess has been used each 500 iterations. The test image size was 256 x 256 pixels. The other 
parameters of the computation were set to 

a = 10 , A = 2, Ax = 1 and At = 10"^ 

and the blur kernel is a Gaussian kernel of standard deviation 5 pixels. 

The computed image is quite satisfying (see figure 7.) However, we note a bad reconstruction 
along the symmetry axis due to the problem geometry and a lack of information. Moreover we 
have to improve the algorithme behavior. Indeed, we observe numerical instability (in spite of the 
regularization process) that leads to a very small time step choice. Therefore the computational 
time is quite long (about 2.5 hours). In addition, classical stopping criteria are not useful here : 
the expected solution corresponds to a "flat" level of function $ and the difference between two 
consecutive iterates means no sense. An estimate of the cost function decrease is not appropriate 
as well (we observe oscillations). We decided to stop after a large enough number of iterations 
(here 20 000). 

In spite of all these disadvantages, this method is satisfactory considering the low signal to noise 
ratio of the radiograph. These good results can be explained by the strong assumptions that we 
add (in particular the binary hypothesis) which are verified by our synthetic object. Anyway, the 
method has beeen successfully tested on "real" images as well, that is imgaes of objets with the 
same kind of properties ("almost" binary) but we cannot report them here (confidential data). A 
semi-impUcit version of the algorithm is actually tested to improve stabiUty. 
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Figure 7. Experimental results 
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